AO-A045  587 


UNCLASSIFIED 


DAVID  « TAYLOR  NAVAL  SHIP  RESEARCH  AND  DEVELOPMENT  CE—ETC  F/6  20/11 
NATFREO#  a computer  PR06RAM  FOR  CALCULATING  THE  NATURAL  FREQUEN — ETC(U) 
JAN  77  J B WILKERSON 

OTNSRDC/ASEO-370  ML 


AO 

A045M7 


IIMMIIII 


NATFREQ  - A COMPUTER  PROGRAM  FOR  CALCULATING  THE 
NATURAL  FREQUENCY  OF  ROTATING 
CANTILEVERED  BEAMS 

by 

Joseph  B.  Wilkerson 


APPROVED  FOR  PUBLIC  RELEASE:  DISTRIBUTION  UNLIMITED 


AVIATION  AND  SURFACE  EFFECTS  DEPARTMENT 


ASED  370 
January  1977 


DAVID 


TAYLDR 

NAVAL 

SHIP 

RESEARCH 

AND 


lllVi»ll!l7Jl»l 


UNCLASSIFIED 


security  classification  of  this  page  r<Win  D«»«  Enitfd) 

REPORT  DOCUMENTATION  PAGE  BEFORE^coMPLETmo^oRii 

l—«rpnNT  niimBE~  / 2.  GOVT  ACCESSION  NO.  1.  RECIPIENT'S  CATALOG  NUMBER 

/7ynTNSRnr7ASED-37flf  I 


'DTNSRDC/ASED-37( 


A TjTLg  fMirf  

( / natfr^/a  computer  program  for  calculating 

\J^  TOE*‘HatWl  EREQUENCY  of  jotating^ 


cantIlevered^eam^ 


l7.  *dTH6Rf«; 


9.  TYPE  OP  REPORT  * PERIOD  COVERED 


9.  PERFORMING  ORG.  REPORT  NUMtCR 


I •.  CONTRACT  OR  GRANT  NUMiCRr*) 


Joseph  Wllkerson 


*.  PERFORMING  organization  NAME  ANO  ADDRESS  10 

David  W.  Taylor  Naval  Ship  R&D  Center 
Aviation  and  Surface  Effects  Department*^ 

Bethesda,  Maryland  20084 

II.  controlling  OFFICE  NAME  AND  ADDRESS  _ ^ 

Naval  Air  Systems  Command  ( 

AIR-320D 

Washington.  D.C.  20361 j 

14.  MONITORING  AGENCY  NAME  9 ADDRESSf'i/  d///«r«nf  from  ControlUng  Otticm)  18 


10.  PROGRAM  ELEMENT,  project.  TASK 
AREA  A WORK  UNIT  NUMBERS 


Program  Element  62211N 
Task  Area  F32.421.210 
Work  Unit  1-1690-100 

IZ.  REPORT  DATE 

Janmaa^F  A77  j 

18.  SBEmlMTM*##WWyT9/  tht»  r0port) 


F3a 


UNCLASSIFIED 

ISa.  OECLASSIFICATION/OOWNGRADINO 
SCHEDULE 


IIS  distribution  statement  fat  ttila  R»por(> 


Approved  for  Public  Release:  Distribution  Unlimited 


D D C 

?f7ar?nn  nr? 

OCT  2T  1977 


I 17.  distribution  statement  fat  lha  tbalrmet  anland  In  Bloc*  30,  II  iHllaranl  tram  Report; 


I IS.  supplementary  notes 


1 1*.  KEY  WORDS  (Conilnua  on  roooroo  old#  II  noeoooarr  and  Idanllly  by  block  nimbtt) 


Helicopter 
Rotor  Blades 


Natural  Frequency 
Structural  Dynamics 


20.  ABSTRACT  (Conilnuo  on  rororoo  oirfo  II  nocoooarp  and  Idantllr  *r  Sloe*  maiAor; 

^ A computer  program  was  developed  to  evaluate  the  natural  frequencies 
of  model  helicopter  rotor  blades  for  use  in  wind  tunnel  evaluations. 

This  program,  NATFREQ,  calculates  the  uncoupled  natural  frequency  and 
corresponding  mode  shape  of  the  first,  second,  and  third  natural  bending 
modes  for  a nonuniform  cantilever  beam  rotating  in  a vacuum.  The  program 
includes  centrifugal  stiffening  effects  and  allows  for  arbitrary  radial  . 

(Continued  on  reverse  side) 


EDITION  OF  I NOV  U It  OBSOLETE 
t/N  OlOa.OIS'CSOI  I 


UNCLASSIFIED 


UNCLASSIFIED 

ul-U*4|TV  CLASSIFICATION  OF  THIS  PAGEOFIlFn  OMa  tnlFratf) 


(Block  20  continued) 


distributions  of  blade  mass  and  stiffness  properties.  The  fundamental 
uncoupled  torsional  natural  frequency  may  also  be  calculated  for  a 
nonrotating  cantilevered  nonuniform  beam.  The  blade  Is  represented  by 
a series  of  concentrated  masses  connected  by  massless  flexures.  The 
method  of  calculation  for  both  bending  and  torsion  Is  by  the  use  of 
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solution  shows  that  good  results  are  obtained  in  the  fundamental  mode 
for  as  few  as  five  mass  elements.  However,  calculated  frequencies  for 
the  second  and  third  modes  of  vibration  show  progressive  error.  This 
suggests  that  up  to  20  elements  should  be  used  for  satisfactory  results 
In  the  higher  modes  for  nonuniform  properties. 
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ABSTRACT 


A computer  program  was  developed  to  evaluate  the  natural 
frequencies  of  model  helicopter  rotor  blades  for  use  In  wind 
tunnel  evaluations.  This  program,  NATFREQ,  calculates  the 
uncoupled  natural  frequency  and  corresponding  mode  shape  of 
the  first,  second,  and  third  natural  bending  modes  for  a 
nonuniform  cantilever  beam  rotating  In  a vacuum.  The  program 
Includes  centrifugal  stiffening  effects  and  allows  for  arbi- 
trary radial  distributions  of  blade  mass  and  stiffness 
properties.  The  fundamental  uncoupled  torsional  natural 
frequency  may  also  be  calculated  for  a nonrotating  canti- 
levered nonuniform  beam.  The  blade  is  represented  by  a series 
of  concentrated  masses  connected  by  massless  flexures.  The 
method  of  calculation  for  both  bending  and  torsion  is  by  the 
use  of  influence  coefficients  and  matrix  algebra. 

Comparison  of  calculated  results  for  a uniform  beam 
with  the  exact  solution  shows  that  good  results  are  obtained 
in  the  fundamental  mode  for  as  few  as  five  mass  elements. 
However,  calculated  frequencies  for  the  second  and  third  modes 
of  vibration  show  progressive  error.  This  suggests  that  up 
to  20  elements  should  be  used  for  satisfactory  results  in  the 
higher  modes  for  nonuniform  properties. 


ADMINISTRATIVE  INFORMATION 

The  work  presented  herein  was  authorized  and  funded  by  the  Naval 
Air  Systems  Command  (AIR-320)  under  Program  62211N,  Task  F32.A21.210 
and  was  accomplished  in  1970.  The  David  W.  Taylor  Naval  Ship  Research 
and  Development  Center  (DTNSRDC)  Work  Unit  was  1-1619-100.  Preparation 
of  this  report  was  funded  under  Work  Unit  1-1619-200. 

All  data  required  for  the  computer  program  described  herein  are 
in  U.S.  customary  units.  Measurements  of  model  geometry  and  frequency 
measurements  were  also  taken  directly  in  U.S.  customary  units.  Hence, 
U.S.  customary  units  are  the  primary  units  in  this  report.  Metric 
units  are  not  given  in  order  to  avoid  confusion  with  the  required 
program  Inputs. 
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INTRODUCTION 

The  design  of  helicopter  rotor  blades  requires  the  ability  to 
analytically  predict  their  natural  frequencies.  It  is  also  necessary 
to  know  the  frequencies  of  model  rotor  blades,  even  If  they  are  not 
dynamically  scaled,  to  avoid  resonant  amplification  conditions.  Resonant 
amplification  occurs  when  the  blade  natural  frequency  coincides  with  an 
integer  multiple  of  the  rotor  rotational  frequency.  The  condition  is  due 
to  the  existence  of  aerodynamic  forcing  functions  which  occur  at  integer 
multiples  of  rotational  frequency.  Operation  at  such  a condition  results 
in  excessive  blade  deflections  taking  the  mode  shape  corresponding  to  the 
frequency  of  excitation.  For  magnitudes  of  forcing  functions  which  are 
low  relative  to  blade  stiffness,  this  results  only  in  higher  stress 
levels;  but  for  higher  magnitudes  of  forcing  functions  (occurring  at  or 
near  the  blades  fundamental  frequency) , the  resonant  amplification 
dominates  the  blade  motion  causing  severe  fatlque  stresses  as  well  as 
affecting  aerodynamic  performance  data. 

COMPUTER  PROGRAM 


GENERAL  DESCRIPTION 

The  computer  program  NATFREQ  was  developed  by  the  author  to  satisfy 
the  above  requirements.  This  program  calculates  the  uncoupled  natural 
frequency  and  corresponding  mode  shape  of  the  first,  second,  and  third 
natural  bending  modes  for  a nonuniform  cantilever  beam  rotating  in  a 
vacuum.  The  program  includes  centrifugal  stiffening  effects  created  by 
spinning  of  the  rotor  blade  and  allows  for  arbitrary  radial  distributions 
of  blade  mass  and  stiffness  properties.  Since  this  analysis  is  intended 
for  rotor  blades,  henceforth  the  beam  length  will  be  referred  to  as 
radius  and  intermediate  lengths  along  the  beam  will  be  referred  to  as 
radial  positions.  In  helicopter  rotor  terminology  the  bending  of  the 
beam  is  taken  to  be  in  the  out-of-plane  direction.  That  is,  the  blade 
bending  takes  place  in  a plane  perpendicular  to  the^  plane  of  rotation 
of  the  rotor  as  shown  in  Figure  1.  Conditions  of  precone  and  aerodynamic 
dynamic  damping  are  not  accounted  for. 
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The  first  uncoupled  torsional  natural  frequency  and  mode  shape  Is 
also  calculated  for  the  cantilevered  nonuniform  beam.  Rotor  rotational 
effects  on  the  blade  torsional  frequency  have  not  been  Included  in  the 
analysis  since  they  generally  have  small  contributions  at  the  relatively 
high  torsional  frequencies  of  model  rotor  blades. 

ANALYSIS 

The  method  of  calculation  for  both  bending  and  torsion  is  by  the  use 
of  Influence  coefficients  and  matrix  algebra  as  outlined  in  Hurty  and 
Rubinstein.^  The  blade  Is  considered  to  consist  of  a series  of  con- 
centrated masses  connected  by  massless  flexures;  see  Figure  2.  Root  end 
conditions  of  the  blade  are  taken  at  the  radius  r^,  which  Is  adjusted  to 
suit  the  particular  geometry  to  be  modeled.  Blade  mass  locations  r^, 

.*.»  from  the  center  of  rotation  are  arbitrary  and  do  not  have  to 
be  evenly  spaced.  Flexures  between  the  masses  may  have  any  distribution 
of  stiffness,  but  an  individual  flexure  is  considered  to  have  uniform 
stiffness  over  Its  length  for  evaluation  of  the  Influence  coefficients. 

The  accuracy  of  results  therefore  depends  on  the  number  of  masses  and 
flexures  used  in  the  model,  as  more  mass  elements  provide  better  repre- 
sentation of  even  a uniform  beam. 

Representation  of  the  blade  for  torsional  calculations  Is  similar  to 
that  for  bending.  Concentrated  mass  polar  moments  of  Inertia  are  connected 
by  torsional  springs,  as  shown  in  Figure  3.  The  radius  r^  for  the  root  end 
condition  and  the  radial  positions  of  the  polar  moments  of  Inertia  are 
taken  to  be  the  same  numerical  values  as  used  for  the  blade  bending 
calculations.  As  previously  mentioned,  the  torsion  mode  Is  uncoupled  from 
the  bending  mode  and  Is  calculated  only  for  nonrotating  conditions. 


CENTRIFUGAL  STIFFENING 

The  equations  of  motion  for  bending  of  the  nonrotating  blade  may  be 
expressed  In  matrix  form  as 

^Hurty,  W.C.  and  M.F.  Rubinstein,  "Dynamics  of  Structures,"  Prentice- 
Hall,  Inc.,  Englewood  Cliffs,  N.J.  (1964)  pp.  110-136. 


3 


Im]  {6}  + [Kl  {6}  - {0} 

where  [K]  « the  square  symmetric  stiffness  matrix 
[ml  “ the  diagonal  mass  matrix 
{6}  * the  column  matrix  of  deflections 
At  each  mass  element  centrifugal  forces  act  to  restore  the  blade  to 
Its  undeflected  position;  see  Figure  4.  The  effect  of  the  restoring 
force  is  compounded  by  the  moment  created  about  position  and  other 

Inboard  elements.  This  restoring  force,  proportional  to  the  square  of 
rotational  speed,  has  a considerable  impact  on  the  blade  out-of-plane 
bending  frequencies.  The  centrifugal  force  acting  at  mass  element  m^^  is 

CF^  = m^  r^ 

This  force  exerts  a restoring  force  component  perpendicular  to  the  blade, 
shown  in  Figure  4,  whl  approximately 

^ = CF^  sin  0^ 


where  sin  0^  = ^i^''  ’^i^  "^i^" 

reduces  to 


For  small  deflections,  6^^  < < r^,  this 


sin  0^  = 


which  gives  the  restoring  force  to  be 


This  restoring  force  is  added  to  the  matrix  equation  of  motion  to  give 

[m]  {6}  + [K]  {5}  + £2^  [m]  {6}  - {0} 

for  the  governing  equation  of  natural  vibration  for  the  rotating  blade. 

Solution  of  the  equation  of  motion  is  accomplished  by  matrix 
algebra  to  give  the  characteristic  eigenvalue  (frequency)  and  eigenvector 
(mode  shape).  The  method  naturally  gives  the  lowest,  or  fundamental, 
frequency  which  will  satisfy  the  equation  for  natural  harmonic  motion. 

The  second  mode  shape  and  frequency  is  a forced  solution,  requiring  it 
to  be  orthogonal  to  the  first  mode.  The  third  mode  shape  and  frequency 
is  then  forced,  requiring  it  to  be  orthogonal  to  both  the  first  and  second 
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modes.  This  is  accomplished  by  the  sweeping  matrices  of  reference  1. 

Since  the  higher  modes  are  based  on  the  fundamental,  they  reflect 
accumulated  error  and  are  therefore  less  accurate  than  the  fundamental. 

INPUT/ OUTPUT 

There  are  two  options  within  the  NATFREQ  program.  The  first  option 
is  the  most  accurate,  allowing  specific  distributions  of  mass  and  stiffness 
properties  to  be  read  in.  This  requires  that  these  distributions  be 
precalculated  from  detailed  engineering  drawings  or  available  from  some 
other  source.  The  second  option  is  used  for  preliminary  estimation  (or 
for  quick  approximations)  of  frequencies  and  mode  shapes  when  less 
information  is  available.  This  option  requires  only  a minimum  of  inputs 
and  uses  programmed  equations  to  internally  estimate  distributions  of  mass 
and  stiffness.  The  procedure  for  calculating  frequency  is  the  same  for 
either  option  once  the  property  distributions  are  obtained. 

Input  format  for  the  first  option  is  shown  in  Table  1.  The  labels 
shown  in  Figure  5 will  help  to  identify  and  understand  the  inputs  of 
Table  1. 

Input  format  for  the  second  option  is  shown  in  Table  2.  Use  of  this 
option  requires  only  blade  root  (RO)  and  tip  (R)  general  geometry.  Mass 
and  stiffness  properties  are  calculated  by  internal  function  statements 
which  are  defined  at  the  beginning  of  the  program.  Those  in  the  current 
version  calculate  radial  distributions  for  an  elliptical  blade  cross  section 
with  the  ellipse  major  axis  in  the  plane  of  rotation.  Thickness  ratio  and 
chord  length  are  considered  to  be  linear  with  radius  from  root  (RO)  to  tip 
(R) . The  type  of  blade  sections  for  which  this  program  was  written  had  a 
hollow  cross-sectional  shape,  shown  in  Figure  6.  First  approximations  to 
natural  frequency  were  obtained  for  this  shape  by  assuming  solid 
elliptical  properties  and  then  subtracting  those  properties  attributed  to 
the  hollow,  or  cut  out,  cross-sectional  area.  Several  inputs  of  this 
option  are  devoted  to  this  purpose  and  are  described  in  Table  2. 

Use  of  the  second  option  for  beams  of  solid  elliptical  cross  section 
is  easily  done  by  setting  all  of  the  cut  out  inputs  to  zero.  Use  of  this 
option  for  beams  with  cross-sectional  geometry  other  than  elliptical  will 


require  changing  the  program  function  statements  to  equations  which 
are  appropriate  for  the  geometry. 

Groups  of  data  may  be  submitted  back  to  back  with  either  program 
option.  The  program  will  read  the  first  data  set,  execute,  return  to 
read  a second  data  set,  and  so  on.  Normal  exit  is  produced  by  putting 
a card  at  the  end  of  the  last  data  set  with  a zero  in  column  5.  (The 
program  exits  when  the  variable  N is  read  in  as  zero.)  A complete  program 
listing  of  NATFREQ  is  included  in  the  appendix. 

A typical  output  for  each  of  the  two  options  is  shown  in  Tables  3 and 
4.  Table  3 (for  the  first  option)  shows  a printout  of  the  input 
distributions  and  the  Influence  coefficient  matrices  for  bending  deflection, 
bending  slope,  and  torsional  deflection.  The  natural  frequency  and  mode 
shape  follow  for  each  of  the  first  three  bending  modes  and  for  the  first 
torsion  mode,  in  that  order.  The  output  frequency  has  the  units  of 
radians  per  second,  and  the  mode  shape  is  normalized  to  1.0  at  the  last 
mass  point  XD(N)  rather  than  at  tip  R. 

Table  4 (for  the  second  option)  shows  the  printout  of  input  blade 
geometry  and  material  properties,  followed  by  the  calculated  distributions 
of  blade  mass  and  stiffness  properties.  The  rest  of  the  output  is  the 
same  as  for  the  first  option  described  above. 

RESULTS 


BENDING 

A lumped  mass  model  for  calculation  of  natural  frequency  allows  for 
better  description  and  greater  accuracy  of  results  for  the  general  case 
of  nonuniform  mass  and  stiffness  distribution.  However,  some  Judgment  is 
required  in  modeling  the  beam  to  be  analyzed  to  ensure  that  a reasonable 
representation  is  provided  to  the  program.  A nonrotating  beam  with  uniform 
mass  and  stiffness  distribution  was  analyzed  as  a check  case  to  evaluate 
program  accuracy  and  the  sensitivity  of  results  to  the  number  of  mass 
elements  used.  The  results  (shown  in  Figure  7)  indicate  that  20  mass 
stations  are  required  to  give  near  exact  results  for  the  fundamental 
frequency  of  the  uniform  beam;  but  even  5 elements  can  give  very  close 


results.  For  this  same  number  of  stations  the  second  and  third  natural 
frequencies  show  progressive  error,  indicating  that  10  stations  or  more 
should  be  used  for  good  accuracy  in  the  higher  modes.  Figure  8 shows  the 
calculated  mode  shape  for  10  mass  elements  and  for  20  mass  elements 
compared  to  the  exact  mode  shape  for  a uniform  beam.  As  shown,  the 
calculated  mode  shapes  are  in  excellent  agreement  even  for  10  mass 
elements  for  the  first  two  modes.  Note  that  these  modes  have  been  adjusted 
to  be  normalized  to  1.0  at  blade  tip  R for  comparison  to  the  exact 
solution. 

The  above  adjustment  is  easily  done  by  assuming  that  the  mode  shape 
is  linear  near  the  blade  tip  (a  good  approximation  for  the  first  three 
mode  shapes  of  cantilevered  beams).  The  calculated  mode  shape  is  adjusted 
by  linear  extrapolation  of  the  two  end  points  to  estimate  the  deflection  at 
tip  R by 

(R-r  ) 

6=6  + (6  - 6 ,)  7 ^ 

R n n n-1  (r-r  ,) 

n n-1 

where  6^  is  the  last  value  from  the  output  mode  shape  (1.0)  and  is 

the  next  to  last  value.  The  new  mode  shape,  is  then  obtained  from 


6' . = 6./6„ 

1 i R 


The  program  has  been  applied  to  the  calculation  of  several  model  rotor 

blades.  One  such  application  was  for  the  blade  with  the  root  and  tip 

cross-sectional  geometry  shown  in  Figure  6.  These  blades  incorporated  8.6 

degrees  of  geometric  twist,  as  depicted  in  the  figure,  and  were  tapered  in 

thickness  ratio.  The  cantilevered  root  end  began  at  r = 0.333  feet  and 

o 

the  blade  tip  was  at  R “ 3.333  feet.  Mass  and  stiffness  properties  were 
calculated  from  the  detailed  engineering  drawings  for  20  stations.  The 
results  from  NATFREQ  for  these  blades  are  shown  in  Table  5 as  compared  to 
the  blades  experimentally  measured  nonrotating  frequencies.  Two  of  these 
blades  were  manufactured;  the  measured  frequencies  from  both  blades  are 
shown.  Table  5 shows  excellent  to  good  agreement  between  the  measured  and 
calculated  natural  frequencies. 
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TORSION 


Torsional  natural  frequency  was  also  checked  against  the  measured 
frequency  of  the  above-mentioned  blades.  As  shown  In  Figure  6 the  model 
blades  had  a hollow  duct  along  the  entire  blade  length  with  an  open  slot 
at  the  rear  edge  of  the  duct.  This  slot  made  the  blades  behave  similarly 
to  an  open,  or  C,  section.  Partial  closure  was  provided  by  discrete 
posts  (near  the  duct  rear  edge)  which  were  spaced  out  along  the  blade 
radius.  Therefore,  effective  values  of  J estimated  along  the  blade  span 
were  used  In  NATFREQ  for  torsional  frequency  calculations.  The  comparison 
of  calculated  to  measured  torsional  frequency  Is  shown  In  Table  S.  Although 
the  comparison  Is  not  as  good  as  for  blade  bending.  It  is  acceptable  con- 
sidering the  effective  J values  which  were  used. 

Torsional  calculations  were  also  made  for  a uniform  beam  as  a check 
case.  Figure  9 shows  these  results  as  they  vary  with  the  number  of 
concentrated  Inertias  used  to  represent  the  beam.  The  fundamental  mode 
shows  only  0.4-percent  error  from  the  exact  uniform  beam  solution  for  a 5 
inertia  element  representation.  This  suggests  that  only  5 to  10  inertia 
stations  may  be  required  for  good  accuracy  when  calculating  the  natural 
frequency  of  most  nonuniform  beams. 

CONCLUDING  REMARKS 


The  computer  program  NATFREQ  has  been  shown  to  be  quite  accurate  for 
the  cases  of  uncoupled  bending  and  torsional  natural  frequencies  of 
uniform  beams.  Calculated  natural  frequencies  for  a model  rotor  blade 
also  showed  good  agreement  with  measured  values.  Care  should  be  exercised, 
however,  to  ensure  that  the  program  is  provided  realistic  physical  Inputs 
when  modeling  any  beam.  For  Instance,  beams  which  have  sudden  changes  In 
stiffness  or  mass  In  some  region  should  have  more  mass  (flexure)  elements 
In  the  model  for  that  region.  Also,  beams  which  have  a point  mass  at  some 
radius  In  addition  to  their  distributed  mass  should  be  modeled  accordingly. 
An  additional  mass  element  should  be  used  at  the  appropriate  radius  in  lieu 
of  distributing  the  point  mass  Into  the  adjacent  beam  elements. 
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PERCENT  ERROR  FROM  EXACT  SOLUTION 


TABLE  1 - NATFREQ  INPUTS  FOR  FIRST  OPTION 


Format 


Variable 


Description 


15  (one  card) 


15  (one  card) 

5E12.5  (up  to  5 cards) 


5E12.5  (up  to  5 cards)  ARAJ* 


Number  of  mass  stations 
(not  to  exceed  21) 

Option  control  (1  for  this  table) 

Array  of  N section  area  moments 
of  inertia,  ft^ 

Array  of  N section  area  polar 
moments  of  inertia,  ft^ 


5E12.5  (up  to  5 cards) 


Array  of  N discrete  mass  weights, 
lb 


5E12.5  (up  to  5 cards) 


Array  of  N discrete  mass  polar 
moments  of  inertia,  Ib/ft^ 


5E12.5  (up  to  5 cards) 


3E10.4,  3F10.2 
(one  card) 


Array  of  N mass  stations  from 
center  of  rotation,  ft 

Youngs  modulas  of  material,  Ib/ft^ 


Shear  modulas  of  material,  Ib/ft 

Density  of  material,  Ib/ft 

Rotor  rotational  frequency,  rad/sec 

Blade  tip  radius  from  center  of 
rotation,  ft 


Blade  root  radius  from  center  of 
rotation  to  root  end  condition,  ft 


These  values  may  be  set  equal  to  1.0  if  torsional  calculations  are  not 
desired;  they  should  not  be  set  to  zero. 

**Set  G to  a large  number  (i.e.,  10^^)  if  torsional  calculations  are  not 
desired. 


TABLE  2 - NATFREQ  INPUTS  FOR  SECOND  OPTION 


i 


J 


Format 


15  (one  card) 

15  (one  card) 
4F10.2  (one  card) 


4F10.2  (one  card) 


5E12.5  (one  card) 


3E10.4,  F10.2 
(one  card) 


Variable 

N 

10 

TCR 

TCT 

CR 

CT 

R 

RO 

ATP 

ART 

ZMIT 

ZMIR 

ZMJT 

ZMJR 

FJ 

E 

G 

DEN 

OM 


Description 

Number  of  mass  stations 
(not  to  exceed  21) 

Option  control  (0  for  this  table) 

Blade  root  thickness  ratio 
(taken  at  RO) 

Blade  tip  thickness  ratio 
(taken  at  R) 

Blade  root  chord  (in-plane 
dimension),  ft 

Blade  tip  chord  (in-plane 
dimension),  ft 

Blade  tip  radius  from  center  of 
rotation,  ft 

Blade  root  radius  from  center  of 
rotation  to  root  end  condition,  ft 

Area  of  cross-sectional  cut  out 
at  blade  tip,  ft^ 

Area  of  cross-sectional  cut  out 
at  blade  root  radius,  ft^ 

Area  moment  of  inertia  of  cut  out 
area  at  blade  tip,  ft^ 

Area  moment  of  Inertia  of  cut  out 
area  at  blade  root  radius,  ft^ 

Polar  area  moment  of  Inertia  of 
cut  out  area  at  blade  tip,  ft^ 

Polar  area  moment  of  inertia  of 
cut  out  area  at  blade  root  radius, 
ft*^ 

Empirical  correction  factor  to 
obtain  torsional  J term  from  the 
calculated  section  polar  area 
moment  of  Inertia 

2 

Youngs  modulas  of  material,  Ib/ft 

2 

Shear  modulas  of  material,  Ib/ft 

3 

Density  of  material,  Ib/ft 

Rotor  rotational  frequency,  rad/sec 
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TABLE  5 - COMPARISON  OF  CALCULATED  AND  MEASURED 
NATURAL  FREQUENCIES 


Mode 

NATFREQ 

Measured 

of 

Calculated 

(rad/sec) 

V ibrat ion 

(rad/sec) 

Blade  1 

Blade  2 

1st  Bending 

121.9 

118.7 

121.3 

2nd  Bending 

634.9 

603.8 

612.0 

3rd  Bending 

1470.0 

1521.8 

1545.7 

1st  Torsion 


919.9 


1097.7 


1030.4 


APPENDIX 


NATFREQ  PROGRAM  LISTING 


C PROGRftM  NATFRE.Q 

c This  program  calculates  the  first  three  natural  sending 

C FREQUENCIES  CF  A NONUNIFCRH  BEAM  ROTATING  IN  A VACUUM. 

C program  by  J.  6.  HILKERSON,  CODE  1619 

C 

C FUNCTION  STATEMENTS 

C elliptical  CROSS  SECTION,  TAPERED  THICKNESS  RATIO  AND  PLANFORM 
C THICKNESS  RATIO  <X» 

TC(X)  = ITCT-TCR) •<X-R0)/RP  ♦ TCR 
C CHORD  (X) 

CTX)  = ICT-CRI*(X-RJ)/RP  ♦ CR 
C DUCT  INTlRNAL  area  (X) 

AKXI  = (ATP-ART»*«X-RD»/RP  ♦ ART 
C DUCT  AREA  moment  OF  INERTIA  <X»,  FLAPPING 
ZMI(X)  = (ZMIT-ZMIRI* (X-R3) /RP  ♦ ZMIR 
C DUCT  POLAR  AREA  MOMENT  OF  INERTIA  <XI,  TORSION 
ZMJ(X)  = (ZMJT-ZMJR) •(X-ROI/RP  ♦ ZMJR 
C BLADE  MASS  (X) 

AMS(X)  = DEN*<3.1h159*TC<X>*C(X)**2/<,.-AI  (X»  » 

C BLADE  AREA  MOMENT  OF  INERTIA  (X»,  FLAPPING 

AMKX)  - 3.1<,159*TC(XI»*3»C(XI»»i»/6<,.  - ZMHXI 
C BLADE  POLAR  AREA  MOMENT  OF  INERTIA(X),  TORSION 

AMJIX)  = 3. 1<»159»TC(X)*C(X)  *»;»•(  1 .♦TC(X)  •*2)/64.  - ZMJ«X) 

C BLADE  MASS  POLAR  INERTIA  IS  CONSIDERED  TO  BE  AREAPOLAR  INE RT I A *DENSI T Y 
c ENO  function  statements 
C 

DIMENSION  XO  ( 22)  ,AIC(  22,  ;.2i  , Bid  22,22  ) ,CIC(  22,22) 

DIMENSION  ARAI ( 22 ) , ARM ( 22 ) , AR A J ( 22 ) 

DIMENSION  ZIU(22,22) ,ZIL(22,22) 

DIMENSION  A I CP <22, 22)  ,BN0N(22,1)  , 0 YN ( 2 2 , 22 ) , OYNS ( 22 , 22 ) , Y (2 2) , ■ 

1 YN(22) ,YN1S(22) ,SWPM(2,22) ,ARP(22)  1 

COMMON  AICP.8NON  ! 

5 continue  ! 

READ  j,N 

IF(N)  163,16d,6  ] 

6 continue  \ 

GG=32.174  j 

READ  3,10  ’ 

IF(IO)  7,7,8 

8 continue 

READ  2CC,  (ARAKI)  ,1=1  ,N) 

READ  2GC, (ARAJ(I) ,1=1 ,N)  : 

READ  200,(ARM<I),I=1,N) 

READ  2bG, (ARP(I) , 1=1, N) 

READ  20C, 1X0(1),  1=1, N) 

READ  9,E,G,0EN,0M  ,R,RO 

9 FORMAT(3E10.4,3Fld.2) 

HB=C.Q 

DO  202  1=1, N 

202  HB=NBf ARM(l)  ! 

200  FORMAK  5E12.5)  I 

PRINT  1030  i 

1030  FORMAT(lHl)  1 

PRINT  1022, DEN, E,G  , 

PRINT  1024,OM 
GO  TO  12 


7 CONTINUE 


R£AO  2.TCR.TCT.CR«CT 

READ  2«R«R0 ,ATP. ART 

READ  200t  ZH1T,2MIR.ZMJT,ZMJR,FJ 

READ  StEtOtOEN.OH 

PRINT  lC20tTCR,TCT,CR,CT.R0tR 

PRINT  1C22«0EN.E,0 

PRINT  1024. ON 

1020  FORMAT  UH1,25X.IBLAi3E  GE0M£TRY$//SX.  iROOT  T/C  = $«  F6. 3 .1  OX  .S  TIP  T/C 
1 =i«Fb. j/7X,$ROOT  CHORO  :S iF8 . 2, 6X , ST  IP  CHORD  sS. F8 . 2/SX , IROOT  RAD 
2IUS  =i,F8.2,5X,»TIP  RADIUS  =1 .F8 .2 .5X///I 
1022  F0RMAH25X.SMATERIAL  PROPERTI  ESI//5  X.  SDENSI  TY  *$,  F12 . 6.  7X,  IE  =I,1P 
3EI3.5.7X.SG  =i.lPEl3.57//» 

1024  FORMAT (25X.SFLIGHT  COND IT IONSS//5X , IOMEGA  :t.F8.2l 

2 FCRMAT(hF1C.2I 

3 FOR.MAT(I5) 

C AVERAGE  AREA  INERTIAS  AND  MASS/UNIT  LENGTH  - INTEGRATION  BY 
C FOUR  POINT  NEMTON-COTES  FORMULA  (SIMPSONS  3/8  RULE» 

A N = N 

RP  = ^.si.o 

H = RP/  (b.»ANI 

ARAKl)  a (AMI(Ru)«’3.*AMI(R0»HI*3.*AMI  (R0Y2.*HI»AMnR0»3.«HII  /8. 
ARAJdi  a (AMJ(R0  M3.«AHJ  (RO  « H I « 3 . • AH  J ( RO  »2  .*H)  »AMJ(R0»3.  •HI)  /8. 

C FJ  IS  AN  EXPERIMENTALLY  OETERMINEO  CORRECTION  FACTOR  TO  OBTAIN  J FROM  THE 
C CALCULATED  SECTION  POLAR  AREA  MOMENT  OF  INERTIA 
ARAJdi  a ARAJ(1)*FJ 
H a RP/(3.*AN) 

DO  10  I a 2,N 

J*2»I-3 

XO  a J*RP/12.*ANI  ♦ RG 

ARAJd)a(  ANJ4XOI»3.*AHJ(XOfH)  *3,  • AM  J ( X0*2  . *H)  ♦AMJ  (XOy3.«HII  /8. 

ARAJII)  a ARAJ(I)«FJ 

10  ARAI(I)  a (AHI(XOI*3.*AMI(XOYHI*3.*AMt(XOt2.*H)YAMI(XO»3.*H))/8.0 
UBau.y 

DO  20  I =1.N 
JaI-1 

XO  a J*RP/AN  ♦ RC 

ARM(I)  a (AMS( XOI t T . • AMS ( XOYN I «3 . • AMS ( X0Y2. •HI ♦AMS(X0f3.«HI I^RP/ 

1 (8.*ANI 

ARPII)  a (AHJ  (X0d3.*ANJ(X0YH)«3.*AMJ(X0»2.*H>vAMJ(X0«3.*Hn*0EN 
1 /(8.*GG) 

20  HB=He*ARM(II 

XOUl  a RP/(2.*AN)  ♦ RO 
H a RP/AN 
DO  25  1=2. N 
25  XOdlaXOd'll  • H 
12  CONTINUE 
PRINT  998 

998  FORMATdHl.lOX.IBLAOE  CHARACTERISTICS  AT  COLLOCATION  POZNTSS/I 
PRINT  1000 

1000  FORMAT(//  13X.SXS.14X.tIXXIX)i.l3X.SNT(XI  I. 18X. tJ I X) t . IIX  .1 IZZ ( X 
II .POLARII 

DO  30  I al,N 

30  PRINT  1001. XDm.ARAim.ARMIII.ARAJ(I). ARPII) 

1001  FORMAT(6X.F6.3.4E20.4I 
PRINT  999. MB 

999  FORMATi/Ziax.lBLACE  HEIGHT  «t.F10.2l 

C OBTAIN  INFLUENCE  COEFFICIENTS  FOR  01 SPLACEMENT .SLOPE • AND  TORSIONAL 
c deflection 

C K DENOTES  LOADED  STATION 
AN  a N 
DO  80  K=1.N 
1*1 

C1*0.0 
C2*0.0 
GO  TO  65 

60  Cl*Bir.(T.KI*F*(ARAI(Ttll-ARAI(n  I tri 


1 


38 


C2=(AIC(ItK)-BIC(ItK)«X0(I))*e*ltR*I(I>ll-«R*I(n  MCZ 
I = I»1 

65  BICII,K>= IX0(K)«X0(I)-X0II)**2/2.C«Cll/tE*ARAI(I) ) 

A1C(1.K)=(X0<KI*XC(I) •X0(l>/2.0-X0(n ••3/6.r«X0(I I*C1«C2I /C€<ARAI( 
III) 

IFCI-KI  63.70t70 
70  L=I*1 

00  75  J=L«N 

BIC( JtKl :BIC(K,K) 

AIC( J,KI  = A1C(K,K>  *SIC (K,K)*(X0(Jt-X0(KI ) 

75  CONTINUE 
80  CONTINUE 

TORSIONAL  FLEXIBILITY  MATRIX,  INFLUENCE  COEFFICIENTS 
00  400  J=1,N 
00  4dC  1=1, N 
IF(l-J)  402,402,404 
402  ZIUCI.JI  = 1.0 
GO  TO  406 
404  ZIU(I.J)  = O.C 
406  2IL(J,II  = ZlUd,  J) 

400  CONTINUE 

: OELX  = LENGTH  OF  FLEXIBLE  CONNECTION  BETWEEN  MASSES 

; OXX  = LENGTH  OF  ELEMENT  HITH  MASS  PROPERTIES 
OELX  = XOdI  - RO 
OXX  = 2.*  (XOm  -ROI 
00  420  J = 1,N 
00  422  I = 1,N 
422  OYNil.J)  = 0.0 

OYNiJ.Ji  = OELX/(ARAJUI*GI 
ARP(J)  : ARP(J)*OXX 

OXX  = 2.*<xo(j»i»-xaij»)  - OXX 

420  OELX  = X0(J41)-X0(JI 

CALL  ZATRIX( OYN,ZIU,OVNS,N,N,NI 
CALL  2ATRIX(ZIL,0YNS,CIC,N,N,Nt 
: PRINT  ROUTINE 

83  PRINT  1»J5 

10U5  F UKMAI (1H1,30X,SI  NFLUENCE  COEFFICIENT  MATR 
II  C E S i///,40X,*OEFLECTION  - AIC«I,J)l/» 

L K=1 
KP=N/iO 

85  KCN=0 
NN1=1 
Nl=l 
N2  = 1C 

86  IF(KP-NNU  87,88,89 

87  N2=N 

84  KON=2 

GO  TO  89 

88  IF(MODIN,lOn  84,84,89 

89  00  93  1=1, N 

GO  TO  190,91,92) ,LK 

90  PRINT  1J07,(AIC(I,J),  J=N1,N2) 

GO  TO  93 

91  PRINT  l4d 7, (BIC(I,J) , J=N1,N2) 

GO  TO  93 

92  PRINT  1007, (CIC(I,JI , J=N1,N2I 

93  CONTINUE 

1007  FORHAT(5X,10E12.5I 
NN1=NN1«1 
N1>N1»10 
N2=N2*10 

IF(KCN-l)  94,94,95 

94  PRINT  1009 

1009  FORNATI/30X,S-  - • - > - CONTINUING  ACROSS,  INCREASING  j - • - > - 

1 -»/» 

GO  TO  8b 


<J  o u o u o u ouou 


9S  LK  3 LK  ♦ 1 

GC  TO  (97«97,99,101l  tLK 
97  PRINT  1011 

1011  FORMAT(1H1,40X.SSLOPE  - B1CII«JIS/I 
GO  TO  85 
99  PRINT  1013 

1013  FORMAT  I1H1.40X»8TORSION  - ClCdtJtS/) 

GO  TO  85 
101  CONTINUE 

CALCULATION  OF  0 MATRIX 

00  207  Isl.N 
207  ARM(I)=ARH<It7GC 
(AM)  MATRIX 
230  00  235  J>1«N 
DO  237  I=1«N 

237  AIC(I.JI=ARM(J)»AIC(ItJ> 

235  CONTINUE 

IF(OH)  240t240.2C5 
(A  Cl  MATRIX 
205  CONST  = 0M**2 
00  210  KsltN 
DO  220  1*1. N 

220  AICP<I.K)=CONST*AIC(I.K) 

210  CONTINUE 

(I  * AC)  MATRIX 
00  225  K=l.N 
AICP(K,K)*1.0« AICP(K.K) 

225  BNON(K«1)=1.0 

CALL  NATINV<N«1«0£TERM,I0> 

IFUO-1)  245.245.227 
227  PRINT  2111 

:ill  FORMAT (2X.|***»N0  MATRIX  INVERSE  •••*3) 

GO  TO  360 
OVN  MATRIX 
24J  00  242  X*1.N 
00  242  J*1.N 
242  CVNd.JI  = AlCd.JI 
GO  TO  253 
245  00  247  1*1. N 
DO  247  J=1.N 
OVNd.JI  = 0.0 
00  250  K*1.N 

2 50  OTNd.J)  = OVN(I.  J)»AICP(I.K)*A1CIK.JI 
247  CONTINUE 

MATRIX  ITERATION  FOR  OMEGA  AND  MODE  SHAPE 

INITIAL  MODE  SHAPE 
253  MS*1 

255  00  260  1*1. N 
268  Y(I)*DyNd.N)/OVN(N.N) 

C ITERATION 

SLAM*C.0 
KOK=0 

262  00  265  1*1. N 
VNd)=0.0 
00  265  J*1.N 

265  yNd)*OYNdtJ)*V(J)«VN(I) 

ALAMOA  * YN(N) 

C NORMALI2E  MODE  SHA  E 

00  270  1*1. N 
270  VNdl*YMII/VN(N) 

C CONVERGENCE  CHECK 

IFIA8S(‘<  AMOA*SLAMI-.06C1*ALAMOA)  305.275.275 
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27S  SLAN  = AtANOA 
00  280  I>1«N 
280  Vll)>VNIl) 

KQK-KQK«t 

IF(KQK-2ai  262«262«28S 
285  PRINT  1G3A.  KQK 

i03«»  FORHATI1M1//////20X,$***  NO  CONVERGENCE  ATSIJ*  ITERATIONS***!) 

305  CONTINUE 

300  FRQ  s l.a/SQRT(ALAMOA) 

PRINT  1035.MS»FRQ  ,KQK 

1035  FORHATI1H1.20X.»HOOE  SHAPES. I 3.5X , SNAT . FREQ  «S«E 12.5,5XS1TERATI0N 
IS  «A.I3/) 

PRINT  10*.0  , (VN(II,I=1,N) 
lOAO  FORHAT<10</5X.F12.0>) 

GO  TO  I 320. 350. 360.360) .NS 
C 

C SECOND  NODE  SHAPE 
C 

320  HS=2 

C SWEEPING  HATRIX 

SWPH(1.1>=0.0 
00  330  J=2.N 

330  SHPHI1.3)=-YN( J)*ARH( J)/(ARHI1I*VNI1)) 

C SAVE  FIRST  NODE  SHAPE  AND  OVN  HATRIX 
00  335  1=1. N 
TN1S(II=TN( I) 

00  335  J=1.N 
335  0YNSII.J)=DVNI1.J) 

C NEW  OYN  HATRIX  FOR  SECOND  NODE 

360  L=HS-1 

00  365  1=1. N 
OYNII.L)  =0.0 
DO  365  J=HS.N 
OYNII.J)=OYNS<I.J) 

K*1 

362  OYNII.J)  = OYN(I.J)»OYNS(I.X) *SWPHIK.J) 

K*K*l 

IFIK-U  362.362.365 
365  CONTINUE 
GO  TO  255 
C 

C THIRD  NODE  SHAPE 
C 

356  HS>3 

SriPNil.2)=0.0 

SNPH<2.1)=0.0 

SWPH(2.2)=0.0 

ONOHl  = ARH(l  )*IYN<2)*VN1S(U-YN1SI2)*YN11)  ) 

0N0H2  = ARH(2>*( VN1S(2I  *VN (1) -VNIS 1 1 ) *YN( 2) ) 

00  355  J*3.N 

SNPH(1.J)= (VN( J)*YN1S(2)-VN1S(J)*VN(2)) *ARHIJ)/0N0H1 
355  SWPHI2. Jl >(YNI J)*VN1SI1)-VN1S< J) * YN( 1 ) ) *ARH( J) /DNOH2 
C NEW  OYN  HATRIX  FOR  THIRD  NODE 
GO  TO  360 
C 

c torsional  frequencies 
C 

C FIRST  node 
360  HS  = 6 

00  650  I * l.N 
00  660  J « l.N 
660  OYNSII.J)  > 0.0 
650  OTNS(l.l)  « ARP(l) 

CALL  2ATRIX(CIC.0VNS*0VN.N.N.N) 

C INITIAL  HOOE  SHAPE 
500  00  501  1=1. N 
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$01  V(l)  3 OVN{I«N)/OVN(N,N> 

KQK  > 0 
C ITERATION 

$02  00  $04  1 - ItN 
TN(1I  • 0.0 
00  $04  J - l.N 

$04  VN(I>  = 0YN(1,JI*V(JI  ♦ VNII) 

KONV  = 1 

C NORMAtIZE  AND  CHECK  CONVERGENCE 
00  $oa  I : 1«N 
VN(I)  : YNID/YNINI 

iF(Aas(YNm-Y<in>.oooi«YN(ii)  $os«$os*$o6 

$06  KONV  S 0 
$00  CONTINUE 

IF(KONVI  310. $10. 370 
$10  DO  512  I = l.N 
$12  YII>  = YNII) 

KQK  * KQK  « 1 
IF(KQK-30)  $02. $02. $14 
$14  PRINT  1034.  KQK 
370  OENOH  s 0.0 
HSS  » HS-3 
00  460  J s l,N 

460  OENOH  = 0YN(N. J)»YNIJI  * OENOH 
FRQ  > SQRT(YN(N)/OENOH) 

PRINT  1C50.HSS.FRQ.KQK 

10$0  FORHATdHl.lOX.srORSION  HOOEI . 13.  SX. INAT.  FREO  «l .£12. $« SX.  tITERA T 
IIONS  =S.I3/t 

PRINT  1C40.(YN(I».  I s l.N) 

GO  TO  5 
160  CONTINUE 
STOP 
END 


r ( 
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ooo  o n o ooo  o o o nonnooo 


SUBROUTINE  NATZNV  <N1 f HI ,OETERN*I 01 

HATRIX  INVERSION  WITH  ACCOHPANVING  SOLUTION  OF  LINEAR  EQUATIONS 
NOVEMBER  1692  S GOOD  OAVIO  TAYLOR  MODEL  BASIN  AH  NATl 

general  form  of  dimension  statement 

DIMENSION  A I , ),BI  • ),INDEX(  «3> 


DIMENSION  A(22«22l «BI22.t) . INOEX(22,3l 
COMMON  A.B 

EQUIVALENCE  (IROU.JROMI,  (ICOLUM«UCOLUH) , lAHAX,  T,  SNAP) 

INITIALIZATION 

M = M1 
N=N1 

10  OETERH=1.0 
IS  00  2(]  J=1,N 
20  INOEX(J,3)  = C 
30  00  550  1=1. N 

SEARCH  FOR  PIVOT  ELEMENT 

60  AMAXsC.O 
65  DO  135  J=1«N 


IFIINOEXI J.3i-i> 

60. 

105. 

60 

DC  100  K=1.N 

IF(IN0EX(K.3>-1> 

80. 

100. 

715 

• 0 IF  ( AMAX  -ABS  (AU.KIM  85.  100.  100 

05  IR0H=J 
90  IC0LUM=K 

AHAX  s AHS  (AIJ.KM 
lUU  CUNTlNt'E 
105  CONTINUE 

INOEXaCOLUM.3l  = INDEX  (ICOLUM.3)  tl 
260  INOEX(I.ll=IROH 
270  IN0EX(I.2)=IC0LUM 

INTERCHANGE  ROMS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

130  IF  IIRON-ICOLUMI  160.  310.  160 
160  0ETERM=-DETERM 
150  DO  200  L^l.N 
160  SHAP=AI IROM.L) 

170  A I1R0M,L)  = AI1C0LUM.LI 
200  A (ICOLUM.L)sSNAP 

IF(HI  310.  31C.  210 
210  DO  250  L=l.  H 
220  SNAPsBdROM.U 
230  B(IROH.U=BIICOLUM.L) 

250  B(ICOLUn.L>=SHAP 

DIVIDE  PIVOT  RON  BY  PIVOT  ELEMENT 

310  PIVOT  =AlICOLUH.ICOLUM) 

OETERM=DETERM*PIVCT 
330  ACICCLUK.ICOLUMI=1.0 
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3<>0  00  35b  L>ltN 

350  aucoluh.li^aiicolun.d/pivot 
355  XF(HI  380,  380,  360 
360  00  370  L>1,M 

370  B(IC0LUH,L)=8IIC0LUN,U /PIVOT 
C 

C REDUCE  NON-PIVOT  RONS 
C 

380  00  550  Llsl,N 

390  IF(Ll-ICOLUN)  400,  550,  400 

400  I :AIL1,IC01UHI 

420  A (L1,ICOLUH)30.0  I 

430  00  450  L<1,N  ’ 

450  AIL1,L>:A(L1,L)-A(IC0LUH,L)*T 
455  IFtHI  550,  550,  460 
460  00  500  L=1,H 

500  BILl,LI=B(Lt,L)-B(ICOLUH,L)*T 
550  CONTINUE 
C 

C INTERCHANGE  COLUMNS 

C 

600  00  710  1*1, N 
610  L=N*1-I 

620  IF  (IN0EX(L,l)-IN0EXfL,2l)  630,  710,  630 
630  JR0M*IN0EXIL,1> 

640  JC0LUH*IN0EXIL,2I 
650  00  705  K*1,N 
663  3NAP-A<R, JROM) 

670  A(K,JR0HI*AIK, JCOLUN) 

700  A(X,JCOLUHI*SHAP 
705  CONTINUE 
710  CONTINUE 

DO  73C  K * 1,N 

IF4INOEXIK,3l  -1)  715,720,715 
715  10  *2 

GO  TO  740 
720  CONTINUE 
730  CONTINUE 
10*1 

740  RETURN 

C LAST  CARO  OF  PROGRAM 

END 
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DTNSROC  ISSUES  THREE  TYPES  OF  REPORTS 

(II  OTNSROC  REPORTS,  A FORMAL  SERIES  PUBLISHING  INFORMATION  OF 
PERMANENT  TECHNICAL  VALUE.  DESIGNATED  BY  A SERIAL  REPORT  NUMBER. 

(2)  DEPARTMENTAL  REPORTS,  A SEMIFORMAL  SERIES,  RECORDING  INFORMA- 
TION OF  A PRELIMINARY  OR  TEMPORARY  NATURE.  OR  OF  LIMITED  INTEREST  OR 
SIGNIFICANCE,  CARRYING  A DEPARTMENTAL  ALPHANUMERIC  IDENTIFICATION. 

(3)  TECHNICAL  MEMORANDA,  AN  INFORMAL  SERIES.  USUALLY  INTERNAL 
WORKING  PAPERS  OR  DIRECT  REPORTS  TO  SPONSORS.  NUMBERED  AS  TM  SERIES 
REPORTS;  NOT  FOR  GENERAL  DISTRIBUTION. 
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